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I. INTRODUCTION 


This thesis is a continuation of research being done at the Naval 
Postgraduate School. An investigation is underway to study the addition 
of aluminum and other metallic particles to solid rocket propellants in 
the hope of achieving greater performance characteristics. Performance 
is sensitive to the particle size distribution throughout the rocket 
motor and nozzle. Since no data exists to predict the distribution in 
this region of the motor, direct observation has been undertaken to 
provide the needed information. 

Holographic techniques have been employed to capture the dynamics of 
the combustion chambers of small rocket motors while firing. These 
techniques are being refined and upgraded. Concurrently, improvements in 
the processing of the holograms to extract the particle size distribu- 
tions are also necessary. 

Previous work done in processing the images has included the use of 
a Quantimet 720 [Ref. 1], and more recently, the use of an IBM PC/AT 
with special software and hardware [Ref. 2]. The latter method has 
proven to be the best currently available. In this method, holograms are 
taken of the rocket combustion chamber while the propellant is burning. 
The holograms are then reconstructed. The reconstructed images are 
stored digitally and later processed to extract the needed particle size 
data. 

A problem arises during the processing of the digital image, how- 


ever. During laser reconstruction of the hologram, speckle noise is 


introduced. This speckle noise can be of the same size as some of the 
smaller particles of interest, and so cause a faulty count in the number 
of features present. 

In an attempt to reduce the speckle, Reference 2 investigated a 
spinning mylar disk, averaging techniques, blurring and lowpass filter- 
ing. Twenty microns was the best resolution that was achieved by these 
means. Some of the loss in resolution may have been due to the optics 
involved; nevertheless, it became necessary to explore other means of 
filtering out the speckle to get the best possible resolution. 

This thesis explores three speckle reduction algorithms suggested by 
work in the field of synthetic aperture radars. The initial step was to 
write computer programs to support each algorithm. Once the programs 
were running, a comparison was done to find which of the three was best 
at reducing speckle while retaining as much resolution as possible. 

This thesis is divided into six chapters. Chapter Two describes the 
process involved in reconstructing a hologram, describes the Itex/PC 
driver software routines available to embed into uSer programs, and 
provides other background material. Chapter Three gives a description of 
speckle and introduces the "speckle index", a figure of merit in deter- 
mining how much speckle is present. The description of each of the 
speckle reduction filters is given in detail in Chapter Four. Chapter 
Five is where the filters are compared with each other. They are ranked 
on speed, speckle reduction based on speckle index, and performance 
based on histogram and visual means. The concluding remarks are con- 


tained within Chapter Six. 


II. HOLOGRAM IMAGE PROCESSING 


An IBM PC/AT is the heart of the entire process of image manipula- 
tion, from digitizing an image to speckle reduction to the counting of 
the features. Installed on the IBM PC/AT are a PC Vision frame-grabber 
board, ImageAction and Itex/PC software (all by Imaging Technology, 
Inc.) the computer monitor, a video monitor and a video cassette 
EecOrder . 

A video image from the VCR can be "grabbed" by the frame-grabber 
under control of the special software. Each of the picture elements 
(pixels) in the 512x480 array comprising a video image is assigned an 
integer gray level from 0 to 255 by the frame-grabber. Level 0 on the 
gray scale is blackest-black. Level 255 is whitest-white, while values 
in between are various shades of gray. Each pixel is uniquely 
addressable and its gray level alterable, thereby allowing the capacity 


Bovachieve digital filtering. 


A. SOFTWARE 
1. ImageAction 
The ImageAction software is a set of menu-driven routines for 
use with a mouse. A totally closed system, it can perform image 
graphics, image analysis, image processing and filtering. Particularly 
useful are the "grab" routine described in the preceding paragraph, and 


the "histogram" routine which outputs the image statistics (standard 


deviation, mean, variance) along with the image's histogram. More de- 
tails are deschibper a mpheterences ss 
2, liber Ee 

The Itex/PC software performs most of the same functions as the 
ImageAction software.[Ref. 3] An important advantage is that the Itex/PC 
Pascal subroutines are able to be called up in Microsoft versions of 
Pascal, C or Fortran programs. This allows the user a great deal of 
flexibility in programming. 

Itex/PC consists of subroutines that may be called from a main 
program to perform a specific function or to return a certain value. For 
example, if the subroutine 

CALL THRESH (lowcut,highcut) 
is used to threshold an image, all the pixel values between the integers 
"lowcut" and "highcut" will be displayed as a new value of 255 (white). 
All other pixels will be displayed as 0 (black). If, as in this thesis, 
"highcut" is set to 255 and "lowcut" varied by the user during run 
time, a binary image results which represents the black particles on a 
white background. 

Sometimes an Itex/PC routine is treated as a function. An 
example is 

De RPIXEESG ye 
Here, the gray level of the pixel located at (x,y) is read and the value 
is assigned to the integer L. The value of L may then be varied by one 
of the filtering algorithms and the new value written back into the 
pixel at w(x. y) susime 


CALL WPIXEL (x,y,L). 


LO 


im ee ee eee eee 





Some specifics of Itex/PC bear mentioning. A standard television 
screen is made up of a 512x512 array of pixels. Of these, only 512x480 
are actually visible. The lowest 32 rows are hidden out of view and may 
contain other information. Itex/PC labels the upper left-hand pixel as 
(0,0), rather than (1,1) as in most conventional arrays. The lower 
right-hand corner of visible pixels is therefore (511,4/9). This method 
of labeling can cause some confusion, and requires more attention to 
detail when writing programs which need to address specific pixels. 

On the monochrome PC Vision frame-grabber board, four look-up 
tables (LUT’s) may be utilized. A LUT serves to transform pixel values 
before they are displayed. The standard LUT is linear, meaning that the 
output values equal the input values. Certain actions, such as 
thresholding, change the LUT values. The actual values of the pixels in 
the frame-grabber board memory remain unchanged by the threshold, but 
are simply altered by the LUT in such a way that a binary image appears 
on the screen. If it is desired to make the threshold permanent and 
actually change the frame-grabber memory to their new values, the MAPLUT 
subroutine must be called. If not, a call to INITIA to initialize LUT's 
returns the screen to its "before-threshold" likeness. 

3. Performance Comparison 

Reference 2 describes a procedure to extract the particle size 
distribution from a reconstructed holographic image. The individual 
particles were numbered, sized in the x and y directions, and sorted 
into bins depending on their size. Non-round particles were included or 
excluded, at the users option. Due to limitations in the software avail- 


able at that time, only one-quarter of the video screen could be 
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processed in a timely manner. A full screen image was "quad-squished", 
thereby losing every other row and every other column of pixels. 
Obviously, resolution was degraded by a considerable amount. 

When Itex/PC and MATLAB (matrix manipulation software) became 
available, the computer program of Reference 2 was upgraded. Prior to 
Itex/PC, it took approximately 40 minutes to process a "quad-squished" 
image. After Itex/PC, a full screen image was processed in approximately 
eight minutes. Taking into account the "quad-squish", this represents a 


twenty-fold decrease in processing time. 


B. HOLOGRAPHY 

As mentioned in the introduction, holographic techniques are used to 
capture the image of a rocket combustion chamber during firing. The 
pulsed ruby laser utilized in the recording process uses a glass 
diffuser in its illumination path. This diffuser is necessary to cut 
down the presence of schlieren interference fringes produced by the 
thermal and density gradients surrounding the burning particles in the 
rocket motor. [Refs. 4 and 5] 

The reconstruction process is shown in Figure 2.1. The hologram is 
reconstructed using a krypton laser and viewed through the microscope 
using 2x, 4x or 10x magnification. The diffuse light from the recording 
interferes with the reference wave of the krypton reconstruction laser. 
This random interference causes speckle to be introduced into the 
image. By attaching a 0.5 lux camera and a video cassette recorder to 


the microscope, the speckle-corrupted image is recorded and preserved 


Te? 


for later use. More detailed descriptions of the reconstruction process 


maeein Reterences 4 and 5. 
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Figure’ 2. | Hologram Reconstruction [Ref. 4] 


Once on VCR tape, the desired image is "grabbed" by either 
ImageAction or Itex/PC software as described earlier. The digitized 
image then may be filtered immediately or stored on floppy disk (SAVEFT 
routine) for later use. The next step is to filter out as much speckle 


as possible without reducing resolution. 


C. HISTOGRAMS 

A histogram of an image consists of a bar graph. A separate bar 
exists at each gray level, with the height of the bar proportional to 
the number of pixels at that gray level. Thus, a histogram gives a 


visual representation of the general darkness or lightness of the image 


PS 


and how widely separated in gray level certain features may be. It can 


suggest where the best level to place a threshold might be. 


D. THRESHOLDING 

Ideally, the feature data in an image would be very different from 
the background in gray level, as in Figure 2.2. A threshold could then 
be placed at a level midway between them with the result that the 
feature particles would be black and the background white. The program 
developed in Reference 2 could then be used to properly size and count 


the particles of interest. 


Number Background 
of pixels 


Grey level 


Pigure az Ideal Histogram 


In practice, Figure 2.3 is much more likely. Particles cover much of 
the range of gray levels, as does the background speckle. A threshold 
placed in an unfiltered image would cause portions of some particles to 


disappear while part of the speckle would appear as particles. 
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Figure 2.3 Histogram Showing an Overlap of Gray Levels 


What is required is a method to process out the unwanted speckle 
noise in the background. With a clear separation between the particles 
of interest and the background, a threshold could then be applied with 


meccess. 
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The purpose of this thesis was to implement filters for the reduc- 
tion of speckle. Therefore, it was natural to attempt to quantify the 


amount of speckle reduction with each iteration of a given filter. 


A. DESCRIPTION Gr Soe chre 

Speckle has the characteristics of random multiplicative noise in 
the sense that the noise level increases with the average gray level of 
a local area [Ref. 6]. The presence of speckle reduces ones ability to 
resolve fine detail. 

The following description may help to give an intuitive idea of what 
speckle is. [Ref. 6] First, imagine a two-dimensional image, such as a 
television screen, laying horizontally on a flat surface. The gray level 
of each individual pixel is represented by a vertical tower, thus making 
a three-dimensional image. The height of each tower is proportional to 
its gray level; darker colors (lower numbers) being short towers and 
lighter colors (higher numbers) being tall towers. Since there are 
512x480 (245,760) individual towers, the impression from a distance is 
of a landscape with smoothly flowing valleys and hills. 

In this imagined landscape, a dark feature would appear as a wide, 
deep valley. The darker the feature, the deeper the valley would be. 


Similarly, a light feature would rise from its surroundings like a 


skyscraper. 
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Now, short, narrow, winding gullies and ridges are added to the 
landscape. Wormy and random in appearance, they corrupt the smoothness 
of the surroundings. This is speckle noise. 

Figure 3.1 shows an image with and without speckle present. It can 
be seen how the speckle is easily confused with the smaller feature 


particles. 


Bee oPECKIER. INDEX 
As a figure of merit in determining the amount of speckle reduction, 
the "speckle index" was used. In Reference 7, Crimmins showed that the 
ratio of local deviation to local mean was a reasonable measure, due to 
the multiplicative nature of speckle noise. 
1. Algorithm 
To compute the speckle index, a fortran subroutine was written. 


The algorithm used was: [Ref. 7] 





Tee ney, 
speckle index = aTN d. 23 ae ems 
where: 
oe 
TE se p(r+a,y+d) C32) 
a,b=-1 
= — b 
dev _ymax_, P(z + a,y + 6) _ min Pletayt ) (es) 


In the above equations, M and N are the dimensions of the array 
in the x and y directions, and p(x,y) is the gray level of an individual 


pixel Tocated at x,y. 


uy 





Figure 3.la Sample Image without Speckle 





bigume 3. 1b Same Image with Speckle Noise Added 
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In simple terms, the local deviation is found by subtracting the 
smallest gray level of a particular pixel and its eight immediate neigh- 
bors from the largest gray level of the same nine pixels. The local mean 
is the average of the nine gray levels. 

2. Implementation 

The complete fortran subroutine to calculate the speckle index 
of an image is included in Appendix A. The subroutine was run with 
various values of M and N to determine a good array size. From Figure 
3.2, it can be seen that a small array can distort the calculation if 
there is a sharp local disturbance. The values M = N = 240 were chosen 
as representative array dimensions for speckle index calculations. This 
size represents about 23% of the visible screen, yet does not take an 


excessive amount of time to calculate. 
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Figure 3.2 Speckle Index as a Function of Sample Size 
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3. Performance 
The IBM PC/AT with the math coprocessor returns the speckle 
index to the calling program in 7/5 seconds. 
Values for speckle index range from over 1.0 (rarely) to a 
theoretical limit of zero. The value of zero could only be achieved if 
the entire image was filtered to the point where every pixel had the 


Same gray level. This is, of course, uninteresting. 
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eee or ECKLE REDUCTION FILTERS 


Three filters were implemented in software and will now be discus- 
sed. All of them call subroutines "speckle" [Appendix A] and "trash" 
[Appendix B]} to calculate the speckle index and threshold the image, 
respectively. 

Two filters depend to a certain degree on the statistics of the 


image. The third filter uses a geometric hulling algorithm. 


eo LHE SIGMA BILTER 
l. Description 

The first speckle reduction filter is based on the standard 
deviation of a Gaussian distribution. By definition, 95.5% of all the 
pixels fall within 2 standard deviations on either side of the mean. Any 
pixels within 2 standard deviations of a given pixel’s gray level are 
included in an averaging scheme, whereas those outside the "2-sigma" 
range are excluded. [Ref. 7] 

Obviously, the standard deviation of the image must be known 
beforehand. ImageAction was used to find standard deviations and histo- 
gram plots very quickly and easily. 

If a particular pixel is considerably different from its 
neighbors, perhaps none of the neighbors will be within the 2-sigma 
range. This would indicate a very sharp feature. To avoid the possibil- 
ity that this sharp feature will not be subject to the averaging process 


at all, a cutoff is established. If the total number of pixels inside 
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the 2-sigma range is less than this minimum cutoff number (2 was 

chosen), the four-neighbor average then replaces the central pixel's 
gray level value. 
2. Algorithm 

The sigma filter program is included [Appendix C]. It uses the 


algorithm as follows: [Ref. 8] 


g{z,y) = 


Pee SEG, 5) 6.5) SR 
oe Ba = 2 ° ae 5 
carers 7= ue pee, it Ss = eles), Sa (4.1) 
ee a ioe 7a =y-2 6(t, j) s=z—-29=y-2 
p(z — 1,y) + ple + Ly) + ple 9) os ple ie e el eee 
is) 
where: p(x,y) = gray level of a pixel in 5x5 local array 
g(x,y) = gray level of the filtered central pixel 
1 if (1 — 20)p(z,y) < plt,7) < (1 + 20)p(z, y) (4.2 
6(z,y) = | ce 
0 otherwise 
g¢ = standard deviation 
3. Discussion 


By altering the value of O , varying degrees of filtering re- 
sult. If @ is increased, the sigma range is increased, and so more 
pixels are included in the average. For this paper, however, the true 
standard deviation of each image was used when it was filtered. 

For illustration of the sigma filter's effectiveness, images of 
the Air Force Resolution Target are presented. Figure 4.1 shows an image 
before and after two iterations of the filter. Visually the speckle has 


been reduced considerably, although the edges have been blurred and the 


resolution has decreased somewhat. 
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Figure4.la Unfiltered Image of the Air Force Resolution 
Bar Target (S.1.— 0.30549) 





Figure 4.1b Same Image Filtered with the Sigma Filter 
(Zelceraetons..o...= 0.10278) 
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Histograms of these two images (Figure 4.2) show some separation 
between features and the background speckle after filtering. A threshold 
at 140 is suggested. 

When this threshold is applied (Figure 4.3), the effectiveness 
of the sigma filter is most dramatically displayed. Now it is much more 
evident that good data can be retrieved from the filtered image, whereas 
the excessive speckle in the unfiltered image would lead to faulty data 


detection. 


B. THE LOCADSSIATISEUGS FIPiEr 
The program to implement this filter is included as Appendix D. 
l. Deser ° 
This filter uses local estimates of the mean and variance in a 
5x5 window about the central pixel in question [Ref. 8]. These local 
Statistics are used to calculate a weight, k. The k value then deter- 
mines where the new gray level of the central pixel will be placed: near 
the original value of the central pixel, near the linear average of all 
pixels in the 5x5 array, or somewhere in-between. 
2. Algorithm 
Throughout this algorithm, p(x,y) represents unfiltered pixels 
and g(x,y) represents the filtered central pixel. To compute the 
filtered gray level of a pixel using the local statistics method, the 


local mean, E(g(x,y)}, and variance, var{g(x,y)}), must first be 


estimated: [Ref. 8] 


E(g(t,y)}% se D> DD lis) (4.3) 
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Figure 4.2a Histogram of the Unfiltered Image 





Figure 4.2b Histogram of the Sigma Filtered Image 
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Figure 4.3a Unfiltered Image Thresholded at 140 





Figure 4.3b Sigma Filtered Image Thresholded at 140 
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var{p(z, y)} + E*{9(z, y)} 


Ze ve 
aes E"{9(z,y)} (4.4) 


var{g(z,y)} * 


where: 
1 z+2 yt2 ; 
var{p(z,y)}~ = D> DL, [elt.z) - B{o(z,y)}] (4.5) 


t=2-—2j=>y—2 


The weight k may then be calculated as: 


fe var{g(z, y)} (4.6) 
o? B?{g(z,y)} + var{g(z, y)} 


Finally, the estimate of g(x,y) is: 


CRONE (et ie le(2,y) — Eto(z,y)) | baie 


From close scrutiny of the algorithm, some of the equations may 
be better understood. Equation 4.3 states that the local mean is the 
average of the 25 pixels in the 5x5 local array. Equation 4.5 is the 
variance of the unfiltered pixels in the 5x5 array; the mean value of 
the local array from Equation 4.3 is subtracted from each of the pixels 
in turn. These differences are squared and then summed over the entire 
array. For more detail on how the remaining equations were developed, 
see Reference 8. 

3. Discussion 

Consider a region of an image which is flat, meaning adjacent 
pixels have about the same gray level. Here, the variance approaches 
zero, and so k approaches zero. The new estimate of g(x,y) is therefore 
close to the mean of the local 5x5 array. Consider now a region of high 


contrast, say the edge of a particle. Here, k approaches one and so 
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g(x,y) is close to the value of the central pixel. The filter therefore 
has a marked tendency toward the retention of high contrast edges. [Ref. 


8 | 
Figure 4.4 shows the Air Force Resolution Target after two 


iterations of the local statistics filter, both before and after thres- 


holding at 140. 


C. THE GEOMET REG} beeen 
The program to implement this filter is included as Appendix E. 
1. Description 

Whereas the previously discussed filters depend to a certain 
degree upon the statistics of the image to be filtered, the geometric 
filter is based on non-linear geometric concepts. 

The algorithm was developed in References 6 and 7. It is a one- 
dimensional routine which is run horizontally, vertically and then in 
the two diagonal directions. It applies in each direction a geometric 
hulling algorithm to the image and then to the image’s complement. 

2. Algorithm :[Ref. 6] 
1) Let a=l, b-O (this sets the values for a horizontal run). 


2) Let c=3 (this is a counter to determine when to change 
directions). 


3) Let d=1l (this controls whether the image or its complement 
is being filtered). 


a? g(x,y) = max [p(z,y), min{p(z — a,y- b) = 1,p(z,y) ts 1}| 
for 1 <7 = 47, 1 see (4.8) 


>) 


p(z, y) = max [g(z, y), min{g(z — a,y — b),g(z,y) + 1,9(z+a,y +5) + 1}] (4.9) 
for l= a <i eee Ni 
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Figure 4.4a iiacoonieeceredswith the Local Statistics Filter 
@eiterations,..o.1.—.0.11464) 





Figure 4.4b The Local Statistics Filtered Image 
Thresholded at 140 
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6) If d=1, let a=-a, b=-b, d=0. Go to step 4). 
If d=0, let d=1. Go to step 7). 


’) g(z, y) = min [p(z, y), max{p(z =U b) ot I poe) 7 1}| 


for 1 <2 =< Mio) ee (4.10) 


8) 


p(z,y) = min [9(z, y), max{g(z — 4,y — 6),9(2, 9) = Igitame, y + 5) — 1} ] oe 
for 1 <2 <= VE ls ye 


9) If d=1, let a=-a, b=-b, d=0. Go to step 7). 
If d=0, let d=l. Go to step 10) 


10) If c=3, let a=0, b=l, c=2. Go to step 4).(this sets up a 

vertical run) 

If c=2, let a=l, b=l, c=#l. Go to step 4). (this sets up a dia- 

gonal run) 

If c=l, let a=l, b=-l,c=0. Go to step 4). (this sets up the 

other diagonal) 

11) If c=0, stop. 

3. Discussion 

The algorithm used exhibits the interesting property that it 
will allow a certain curvature in the terrain (referring again to the 
"landscape" mentioned earlier), but will tear down anything more dras- 
tic. With a few exceptions [Ref. 6], a curve sharper than 45 degrees at 
any vertex will be filtered. What this means is that the narrow valleys 
comprising speckle will be filled in and the walls torn down after 
sufficient iterations of the filter. Of course, the objects of interest 
that should be preserved are also filled in and torn down, but at a much 
slower rate. Generally, the larger and darker the feature, the more 
slowly it will be degraded. 


If an object of interest is almost as narrow as the speckle, 


they both will be reduced at about the same rate. Hopefully, the object 
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is much darker than the speckle, so the speckle will be beaten down 
after a number of iterations leaving the object of interest largely 
imeact . 

Figure 4.5a shows the Air Force Resolution Chart after three 
iterations of the geometric filter. Figure 4.5b shows the same filtered 


image thresholded at a gray level of 140. 
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Figure 4.5a Image Filtered with the Geometric Filter 
(3 lteratrons. “Sia 0712013) 





Figure 4.5b Geometric Filtered Image Thresholded at 140 
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V. COMPARISON OF THE FILTERS 


This chapter will compare the three filters using the speckle index 
as a figure of merit in gauging how much speckle has been removed. The 
filters will be compared in how much speckle they remove per iteration, 
the time per iteration, histogram differences, and, of course, visual 


@iftferences. 


A. SPECKLE INDEX COMPARISON 

All three filters reduce speckle at their own rate. Figure 5.1 shows 
this fact. The sigma and local statistics filters both reduce speckle by 
almost a factor of two after only one iteration, but then dramatically 
"slow down" with increased iterations. This may be attributed to the 
fact that both filters are actually various themes on the blurring 
technique. More iterations serve to merely increase the blur. 

Figure 5.1 shows that the geometric filter’s curve is more gradual, 
allowing the user to stop at a less severe level of filtering. This is a 
particularly desirable trait, because resolution degradation becomes a 
problem after only a few iterations. It should be noted that filtering 
cannot stop at an arbitrary value of speckle index, only at discrete 
levels depending on how many iterations have occurred. It is foreseeable 
that a given filter may allow a degree of speckle reduction not 


achievable by the others. 
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Figure 5.1 Speckle Reduction Comparison 


On the basis of Figure 5.1, the geometric filter ranks first in its 
ability to control the amount of speckle reduction, followed by the 


local statisti@smer) tere 


B. TIME COMPARISON 


The times for an iteration of each filter are in Table 4.1. 


TABLE 4.1 TIMES PER ITERATION 


FILLER TIME 
Sigma 7 min, 20 sec 
Geometric 18 “mim, S0esec 
Local SGatistices 39 min, 745 sec 
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The local statistics filter is slow due to the nested do-loops, 
frequent calls to subroutines and large number of calculations required 
inside each do-loop. 

It. should be noted that speed is of secondary importance since 
there is no requirement in this application to achieve results within a 
given time limit. Also, dedicated image-processing equipment can do the 
calculations much more efficiently and faster. For example, Reference 8 
reports times of 14 seconds per iteration for the geometric filter using 


a VAX 11/780 and DeAnza IP-5500 digital video processor. 


C. VISUAL COMPARISON 

The ultimate test of any procedure dealing with images is how the 
image looks to the viewer. Comparisons of this type are difficult to 
justify, for each person sees an image slightly differently. Neverthe- 
less, based on numerous runs of the filters, the geometric filter is 
judged best at retaining the edges, basic shape and size of objects of 
interest while beating down the speckle. This can be seen by close 
comparison of Figures 4.lb, 4.4a and 4.5a. All three Figures are close 
in speckle index, yet the local statistics and sigma filters tend to 
smear and spread particle edges more than does the geometric filter. 
This observation has been verified in runs using other images as well. 

Once the threshold at 140 is applied (Figures 4.3b, 4.4b and 4.5b), 
the sigma filter is definitely seen as the inferior of the three. The 
other two are very close. Examination reveals a few small particles of 
speckle that the geometric filter failed to eliminate which the local 


statistics filter successfully handled. This is a consequence of the 
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filtering stopping at discrete levels. In most other cases, the 


geometric filter is superior. 


D. HISTOGRAM COMPARISON 

All three filters produced histograms generally similar in shape and 
distribution (Figures 4.2b, 5.2a, 5.2b). Notable differences include the 
geometric filter’s reduction of the spike at 255, the way it retained 
virtually intact the darker features and its more pronounced "valley" 
between the two "humps". 

The claim was made earlier that the geometric filter had a tendency 
to reduce smaller (hopefully speckle) particles much faster than the 
dark features. The retention of the darker pixels in the histogram is 


evidence that this claim is valid. 
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Figure 5.2a Histogram of Local Statistics Filtered Image 





Figuré«5.2b Histogram of Geometric Filtered Image 
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All the filtered histograms pointed to a threshold at about 140. 
Thus, this was the value used on all the images, including the original, 


for comparison purposes. 
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VI. GONCLUSIONS 


The IBM PC/AT with dedicated software and hardware is a viable 
system to reduce speckle in reconstructed holograms. The ITEX/PC soft- 
ware with the PC Vision frame-grabber board may be used in conjunction 
with fortran programs to achieve the filtering desired. 

While all three filters are superior to current techniques, the 
geometric filter has been found to be the best of the three. It combines 
the ability to hold the edges and shape of objects with the desirable 
trait of less filtering per iteration. In this way, the user has more 
control over the amount of filtering to be done and to what degree the 
speckle will be reduced. 

While the geometric filter has been found to be the best overall 
filter, it may not be superior in every particular circumstance. It has 
been noted that due to the discrete jumps in speckle index which occur 
after every iteration, one filter may be able to reach the region of 
optimum filtering when another cannot. For this reason, the local 
Statistics filter should be used along with the geometric filter and the 
results after 2 or 3 iterations compared. 

Resolution of the original image was about 12 microns. This was 
degraded to about 14 microns in the filtered images. It is doubtful that 
much more improvement can be made toward reducing the speckle degrada- 
tion due to filtering. Any further resolution improvements most likely 
will have to occur in the recording and reconstruction process, and 


better optics. Some improvements have been made recently in this area. 
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Reference 4 reports resolutions typically of 8 microns, and in certain 
cases, 4 microns. Taking into account the degradation after filtering, a 
conservative estimate of 10 microns resolution is now possible in the 


finished filtered and thresholded image. 
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APPENDIX A. 


SPECKLE INDEX SUBROUTINE 


Sinclude: ‘itexpce.inc’ 


OOO OO OOO OOO "On O° AO A A 0 


Q 


welts 
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subroutine speckle(spklindx) 

KEKE HHH HEHE AHH HHH HANH HHH AK AHHH HHH HHH HN HHH HE EE EE ee Seah ae a ek oe ae ak ak 
This subroutine computes the speckle index (a figure of 

merit) of the upper-left hand corner of the image, a 240x240 


array. 
SEK HHH HHH KH HK HH HH HER HH HRI IE EIR IE IE IE IEE aE sea re ak ae de i ae aba sab sia ak ak ae ab ak sea 


**k*X DEFINITIONS *** 


small min gray value in 3x3 local array 

big max gray value in 3x3 local array 

sum summation of the nine gray values in 3x3 array 
dev difference between small and big 

mean average of the nine gray values 

local local deviation divided by the local mean 
total summation of local over entire input array 


spklindx speckle index 
**&* DECLARATIONS *** 


implicit integer*2 (a-z) 
integer*2 small,sum,big,dev 
real*4 mean,local,total,spklindx 


total=0.0 
do 40 m=1,238 
do 30 n=1,238 


*** COMPUTE THE LOCAL MEAN AND DEVIATION *** 


big=max0(RPIXEL(m-1,n-1) ,RPIXEL(m-1,n),RPIXEL(m-1,n+1), 
RPIXEL(m,n-1),RPIXEL(m,n) ,RPIXEL(m, n+l) ,RPIXEL(m+1,n-1), 
RPIXEL(m+1,n),RPIXEL(m+1,n+1) ) 
small=minO(RPIXEL(m-1,n-1),RPIXEL(m-1,n),RPIXEL(m-1,n+1), 
RPIXEL(m,n-1),RPIXEL(m,n) ,RPIXEL(m,n+1),RPIXEL(mt+l1,n-1), 
RPIXEL(m+1,n) ,RPIXEL(m+l1 , n+1) ) 
dev=big-small 
sum=RPIXEL(m-1,n-1)+RPIXEL(m-1,n)+RPIXEL(m-1,n+1l1)+ 
RPIXEL(m, n-1)+RPIXEL(m,n)+RPIXEL(m, n+1)+ 
RPIXEL(m+1,n-1)+RPIXEL(m+1 ,n)+RPIXEL(m+1,n+1) 
mean=sum / 9.0 
if (mean .eq. 0) then 
mean=1 

end if 
local=dev / mean 
total=total + local 

continue 


41 


40 


continue 

spklindx=total/56644. 

write(*,*)’the speckle index is:',spklindx 
return 

end 
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APPENDIX B. 
THRESHOLD SUBROUTINE 
Sinclude:’itexpe.inc’ 


Pub momeime tbasnm 
So tk oink oping ok ok sk oink ob bk ok ok oioinksk ok ok ak akin ok ok ok sk oko ok ok ok ok ok seak ok ok ak ok oko sk ok sk sk skse sles okt 


Cc 
c This subroutine thresholds the current image to whatever 
a value the user decides. It will change this threshold value as 
eC many times as is necessary. 
Cc AK HK HAHAHA HAH HHA AHH HK IER HWE HE TERM WHEN MH HEH METER REDE ae ge at ak ak sea 
Cc 
@ *** DEFINITIONS *** 
c errval error value for the file retrieval 
€ inet threshold level 
‘o 
C ***x DECLARATIONS **%* 

implicit integer*2 (a-z) 

integer*2 errval,limit 
c 
© **x*k SET BASE ADDRESS OF HARDWARE REGISTERS, INDICATE A PSEUDO- 


a COLOR BOARD AND INITIALIZE THE LUT's TO STANDARD VALUES. *** 
call SETHDW(16#f£00, 16#c000,1) 
Calle Initia 


c 
20 write(*,*) 'Thresholding-- enter the desired pixel value’ 

read(*,*) limit 

C 

@ *%*%* THRESHOLD THE 4 LOOK UP TABLES (LUT'S) *** 
Call SELLULGORO) 
Gall THRESHA imi, 255) 
Call TSELLUECIH UO) 
Callies LHR ESC imine 255) 
eae SELLUERCZ 0) 
Galli THRESH Climnes:2 oo) 
Call SELLUT (3-0) 
Cale THREOHC Limit, 255) 

C 

o *** QUERY USER ON WISHES TO THRESHOLD AGAIN *** 
write(*,*) ‘do you wish to set a new threshold value? (0 for yes)’ 
read(*,*) d 
if (d .eq. 0) then 

call INITIA 
go to 20 

end if 

@ *** TNITIALIZE THE IMAGE BEFORE RETURNING **%* 
cal iNt EEA 
return 
end 
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APPENDIX C. 
TWO-SIGMA FILTER PROGRAM 


Sinclude:’itexpc.inc’ 


SRI III III I III II III Inn Eni ini ik ark bik sea sea ok ae aba ak ak 
This is a statistical filter known as a "2-sigma". It calcu- 
lates the range of gray levels within 2 standard deviations of 
the central pixel in a 5x5 array. Any of the 24 pixels surround- 
ing the central pixel that fall within this "2-sigma" range are 
then included in an averaging algorithm. If there are no pixels 
within the given ranges, then the 4 neighbor average is calcu- 


lated to correct for this "sharp-shot noise". 
TORII II III III III III III III Ee ir iki kb bok b ak kkk ek 


*** DEFINITIONS *** 


fname filename 

comlin comments line for file storage 

dev standard deviation 

devo dev converted for program’s use 

less 2 standard dev's below 1, expressed as % 
more 2 standard dev’s above 1, expressed as % 
flag filter iteration number 

delta indicates if a given pixel is within 2dev's 
zee summation of all gray levels within 2 dev's 
E gray value of central pixel in 5x5 array 
small lower 2-sigma limit 

large upper 2-sigma limit 

up pixel directly above central pixel 

dn pixel directly below central pixel 

g(x,y) filtered pixel 

spklindx speckle index 

trash thresholding subroutine 

speckle speckle index subroutine 


%ke& DECLARATIONS *** 


Oo O° O-O8°): O° Oo OA 8 0 8 2 0 0.0 OG 6. eG 0.60 QO 2 f 2 2 2 eae.) 


implicit integer*2 (a-z) 

integer*2 delta,errval,small,large,up,dn,zee,yy,xx, flag 
integer*2 q(5) 

real*4 dev,devo,less,more,spklindx 

character*1l2/7 comlin 

character*20 fname 

integer*2 p(512,240) 


‘e) 


**k* SET HARDWARE REGISTERS, INDICATE USE OF A PSEUDOCOLOR BOARD 
Cc AND INITIALIZE LUT’s TO STANDARD VALUES *** 

call SETHDW(16#FFOO, 16#C000, 1) 

call INITIA 
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Q 


Or OOO © 


Q 


5401 


500 


10 
20 


write(*,*) ‘enter the name of the file to retrieve:’ 
read(*,5401) fname 
format(a) 
errval=READFT(0,0,512,512, fname, comlin) 
write(*,*) ‘error code from READFT=’ ,errval 
if (errval .eq. 0) then 
write(*,*) comlin 
endif 


call TRASH 

write(*,*) ‘enter the standard deviation of the image: ’ 
read(*,*) dev 

devo=dev/256. 

*** COMPUTE SPECKLE INDEX *** 

call SPECKLE(spklindx) 


less = 1-2*devo 
more 1+2*devo 


KKKK KK eke KK Rokk ekki k bk kb ook et 


This is the actual sigma filter 
FAA HHA AA HHA NEHER HER AR HAR HEAR ER HK TE BE 


0 
flag+l 
do 40 y=2,240 
do 30 x=2,509 
delta = 0 
zee=0 
r=RPIXEL(x, y) 
small=nint(less*r) 
large=nint(more%r ) 
do 20 k=y-2,y+2 
call RHLINE(x-2,k,5,q) 
do 110 lel) 5 
Tree) ce. small) sand. (q(1) 
large)) then 
delta=delta + l 
zee = zee + q(1) 


flag 
flag 


ence lr 
continue 
continue 


wee SHARP SHOT NOISE IS CORRECTED *** 
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ves 


30 
40 


OQ 


50 
60 


OQ 


70 
80 


90 
100 


if (delta .le. 2) then 
up=RPIXEL(x, y-1) 
dn=RPIXEL(x, y+1) 


g(x,y)=(uptdntq(2)+q(3)+q(4))/5 


go to 30 
end if 
g(x,y)=zee/delta 


continue 
continue 


**%* WRITE FIRST HALF INTO FRAME-GRABBER *** 


do 60 y=2,240 
do 50 x=2,509 


call WPIXEL(x,y,g(x,y)) 


continue 
continue 
write(*,*)’ 


the first half is frleered” 


*x* FILTER THE SECOND HALF *** 


do 100 y=240,480 
do 90 x=2,509 


delta = O 
zee=0 
r=RPIXEL(x, y) 
small=nint(less*r) 
large=nint(more*r) 
do 80 k=y-2,y+2 

call RHLINE(x-2,k,5,q) 

do 70 1=1,5 


if ((q(l) .ge. small) .and. 


large)) then 
delta=delta + 1 
zee = zee + q(1) 
end if 
continue 
continue 


**x* SHARP SHOT NOISE IS CORRECTED *** 


yy=y -240 

if. (delta: .le- 7) then 
up=RPIXEL(x,y-1) 
dn=RPIXEL(x, y+1) 


g(x, yy)=(uptdntq (2)+q(3)+q(4))/5 


go to 90 
end if 
g(x, yy)=zee/delta 


continue 
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(q(1) 


le. 


110 
120 


0 


400 


*x*%* WRITE SECOND HALF BACK INTO FRAME-GRABBER *** 


do 120 y=240,480 
yy=y-240 
do 110 x=2,509 
eee LIE Ai eal GE; CaN roe eal) 
continue 
continue 
write(*,*) ‘the second half is finished’ 


call SPECKLE(spklindx) 
call TRASH 


**x*k ADMINISTRATIVE *** 


write(*,*) flag,’ =flag’ 

write(*,*) ‘do you wish to store the image? (0 for yes)’ 
read(*,*) o 

if (o .ne. 0) goto 400 

write(*,*) ‘enter filename you wish to store the image under’ 
read(*,5401) fname 

write(*,*) ‘enter any comment, terminated with a !:’ 
read(*,5401) comlin 
errval=SAVEFT(0,0,512,512,1,fname,comlin) 

write(*,*) ‘error code from SAVEFT=’ ,errval 

ieee ban LemeO ye coto. 500 

stop 

end 
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APPENDIX D. 
LOCAL STATISTICS FILTER PROGRAM 


Sinclude: “itexpe. ine. 
KEK KKH HER KKK ETE TERE HRI IE IEEE TE ETE EEE IIE TEE TEE IIE TE DE DES EI ae ae ae ae ai ak aie ai asia 


Cc 

re This is the local statistics filter. It computes the local 
c variance and mean and then uses these numbers to calculate an 
Cc estimate of the pixel value. 

C SKK HATER HAH KKK KK HARK HAR HHA KKH KHER KHER KEE KHER EEK OEE EE 
Cc 

¢ *k* DEFINITIONS *** 

Cc fname filename 

Cc comlin comments line for file retrieval 

ec sum summation of gray levels for 5x5 local array 

eC errval error value for use with file storage 

c flag filter iteration number 

re q array to hold five pixel values 

Cc meanz average of 5x5 local array 

e meanz2 meanz squared 

c sum2 summation of variances of local array 

Cc dev standard deviation, entered by user 

Cc stddev2 dev squared 

Cc Varx variance of image without noise 

€ varz variance of image with noise 

c k2 k parameter 

Cc g(x,y) filtered pixel gray values 

Cc trash threshold subroutine 

es speckle speckle index subroutine 

Cc spklindx speckle index 

Cc 

c *** DECLARATIONS *** 

Cc 


implicit integer*2 (a-z) 

integer*2 sum,errval,flag,q(5) 

real*4 meanz,sum2,meanz2,varz,varx,stddev2 ,dev,k2 
character*127 comlin 

character*21 fname 

integer*2 g(512,402) 


QD 


***k SET HARDWARE REGISTERS, INDICATE USE OF A PSEUDOCOLOR BOARD 
c AND INITIALIZE LUT’s TO STANDARD VALUES x*** 

call SETHDW(16#ff£00 , 16#c000,1) 

call INITIA 


write(*,*) ‘enter the name of the file you wish to call up:’ 
read(*,5401) fname 

format(a) 

errval=READFT(0,0,512,512,fname,comlin) 

write(*,*) errval,’=error code from readft’ 
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QD 


7) 


O? 


‘@) 


0 
Z0 


30 
40 


50 


i: Cereval seq, 0) then 
WeiteC<)*) comlin 

end if 

call TRASH 


write(*,*) ‘enter the standard deviation of the image: ’ 
read(*,*) dev 

stddev2=(dev/256.)**2 

write(*,*) stddev2,'=stddev2' 

flag=0 


*** COMPUTE THE SPECKLE INDEX *** 
call SPECKLE(spklindx) 


flag=flagt+l 
sum=0 . 
sum2=0 . 
do 60 y=2 ,400 
do 50 x=2,509 


*** COMPUTE THE LOCAL MEAN *** 


do 20 l=y-2,y+t2 
cal RHEENE(x-2,1,5,q) 
do 10 k=1,5 

sum=sumt+q (k) 

continue 

continue 

meanz=sum/25. 

sum=0. 


*** COMPUTE THE LOCAL VARIANCE *** 


do 40 l=y-2,y+2) 
call RHLINE(x-2,1,5,q) 


do 30 k=1,5 
sum2=sum2+(q(k) -meanz) **2 
continue 
continue 


TAs eCOMPUDE THE ESTIMATE OF pcx ,Y) *** 


varz=sum2/25. 

sum2=0. 

meanz2=meanz**2 
varx=abs((varz+meanz2)/(stddev2+1) -meanz2) 
k2=varx/(meanz2*stddev2+varx) 

g(X,y)=int (meanz+k2* (RPIXEL(x, y) -meanz) ) 


continue 
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QO 


60 


70 
80 


100 


continue 
*** WRITE g(X,Y) ARRAY INTO FRAME-GRABBER *** 


do 80 y=2,400 
do 70 x=2,509 
call WPIXEL(x,y,g(x,y)) 
continue 
continue 
call SPECKLE(spklindx) 
call TRASH 


write(*,*) ‘do you wish to save to a file? (0 for yes)’ 
read(*,*) O 
if (O .NE. @) yeeros 


write(*,*) ‘enter the filename to store the image under: ’ 


read(*,5401) fname 

write(*,*) ‘enter any comment, terminated with a !:' 
read(*,5401) comlin 
errval=SAVEFT(0,0,512,512,1, fname ,comlin) 

write(*,*) errval, ‘'=error code from SAVEFT’ 

if ({tlag .le, 10) egocons 


stop 
end 
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APPEND IAS: 


GEOMETRIC FILTER PROGRAM 


pmc lude: ’ 1ltexpe. ine’ 


eG awe 3 Or 30.0 8) OO. same Oo Oo. O° OO AO f° O 


Q 


5401 


BKK HHH HHA HAHAH HHH HA AHH HHH HA HHH HAA KHIR KK AIHA TEE a a ae ae ae a a ae 

This is a program to filter an image through the use of a 
geometric algorithm. It calls subroutine ‘speckle’ to find the 
speckle index, and subroutine trash to threshold the image. The 


ITEXPC software is used extensively. 
KKK H HHA A AHA AA HHA A HHA AHHH AHHH AHH HAKKAR HHH EE TE IE re a ae re ae a ae eat 


**x* DEFINITIONS *** 


fname filename 

comlin comments line 

flag iteration number of filter 

t a counter; determines which quarter of image 

Ww row number 

a,b used to control direction within algorithm 

Cc used to determine direction 

d controls whether image or complement is filtered 
Spiel ind speckle index 

q array to hold a complete row of pixel gray levels 
errval error value for use in file storage 

pie) image pixel gray levels 

g(x,y) image complement pixel gray levels 

trash threshold subroutine 

speckle speckle index subroutine 


*k& DECLARATIONS *** 


implicit integer*2 (a-z) 
integer*2 errval,flag,q(512) 
real*4 spklindx 

character*127 comlin 
character*20 fname 

iiecevers2 p(l20',512) ,¢(120,512) 


*** SET HARDWARE REGISTERS, INDICATE USE OF PSEUDOCOLOR BOARD 
AND INITIALIZE LUT’s TO STANDARD VALUES *** 

call SETHDW(16#ff£00, 16#c000, 1) 

call INITIA 


write(*,*) ‘Enter the name of the file to retrieve: '’ 
read(*,5401) fname 

format (a) 

errval=READFT(0,0,512,512,fname,comlin) 

write(*,*) ‘error code from READFT = ’ ,errval 


5 


O 


O 


O 


500 


if (errval .,eq. 0) then 
write(*,*) comlin 

end if 

call TRASH 

flag=0 

w=0 

e=0 

write(*,*)flag,’is the run number’ 

flag=flag+l 


*** COMPUTE THE SPECKLE INDEX *** 


call SPECKLE(spklindx) 
do 1500 t=1,4 


**%* CONVERT THE INPUT IMAGE INTO A 512x120 ARRAY, P(X,Y) *** 


do 2 y=1,120 
Call RALINECO we old) 
w=wtl 
do 1 x=1,512 
p(y, x)=q(x-1) 
continue 
continue 


**x*x THIS IS THE START OF THE GEOMETRIC FILTER *** 


can nk 
Co Oar 


d=] 
write(*,*)’working. .. .’ 
do 120 y=1,120 
do 110 x=1,512 
gcy,x)=max0(p(y,x) ,minO(p(y-a,x-b)-l,pCy,x)+1l)) 
continue 
continue 
do 140 y=1,120 
do 130 x=1,512 
pcy,x)=maxO0(g(y,x) ,minO(g(y-a,x-b),g(y,x)+tl, 
g(y+ta,x+b)+1)) 
continue 
continue 
if (d .eq. 1) then 
=-4 
b=-b 
d=0 
go to 102 
end if 


By 


©) 


105 


899 


300 


BOL 


inca eq. 0) then 
d=] 
end if 
do 160 y=1,120 
dow la0 xe? 512 
g(y,x)=minO(p(y,x) ,maxO(p(y-a,x-b)+l1,p(y,x)-1)) 
continue 
continue 
do 180 y=1,120 
do 170 x=l1 ,512 
p(y,x)=minO(g(y,x) ,maxO(g(y-a,x-b) ,g(y,x)-l, 
g(yta,xtb)-1)) 


continue 
continue 
a1 (cd .eqeae. then 
a=-a 
b=-b 
d=0 
£0 co 105 
end if 
Pin deed. 0) then 
=] 
end if 
if (c .eq. 3) then 
a=0 
b=1 
c=2 
Fo toe lo2 
end if 
ieee ede) then 
=] 
b=1 
c=1 
gon con LOZ 
end if 
tf (o..eq. 1) then 
a=1 
b=-1 
c=0 
go to 102 
end if 


*x*x WRITE BACK INTO FRAME-GRABBER AND DISPLAY «** 


donoUleyv—i, 120 
do 300 x=1,512 
q(x-1)=p(y;x) 
continue 
call WHLINE(0,e,512,q) 
e=e+l] 
continue 


a3 


1500 


LoL 


a27, 


continue 
xe ADMINISTRATIVE *** 


call TRASH 

write(*,*)’'Do you wish to save image to a file? (0 for yes):' 
read(*,*) o 

if (o .ne. ©) goto 1501 

write(*,*)’enter the filename you wish to store the image under: ’ 
read(*,5401) fname 

write(*,*)’enter the comment, terminated with a "!" :' 
read(*,5401) comlin 
errval=SAVEFT(0,0,512,512,1,fname,comlin) 

write(*,*)’error code from SAVEFT = ',errval 

if (flag . le il2)eeote, 500 

call SPECKLE(Cspelindx) 

stop 

end 
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